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This paper develops an algorithmic solution strategy which enables 
handling the positive/indefinite stiffness characteristics associated with 
the pre and postbuckling of structures subject to complex thermomechanical 
loading fields. The flexibility of the procedure is such that it can be 
applied to both finite difference and element type simulations. Due to 
the generality of the algorithmic approach developed, both kinematic and 
thermal /mechanical type material nonlinearity including inelastic effects 
can be treated. This includes the possibility of handling completely 
general thermomechanical boundary conditions. To demonstrate the scheme, 
the results of several benchmark problems is presented. 

INTRODUCTION 

Literally a multitude of studies have been reported on the isothermal 
simulation of problems wherein kinematic and/or material nonlinearity is 
excited. In recent years, most typically such work involved the use of the 
powerful finite element ( FE) scheme [1]. In contrast, much less work is. 
available for nonisothermal versions of such problems. This is an outgrowth 
of several main factors namely: 

i) Unlike mechanical type loads which are generally applied at- specific 
points around a given structure, transient thermally induced loads 
occur at every body point causing complex distributed loading and un- 
loading fields which typically induce difficulties in simulating proper 
inelastic type behavior; 

ii) Since thermal loads are internally induced, for nonlinear situations, 
it is typically quite difficult to adequately forecast the level of 
incrementation necessary for nonlinear equation solvers to yield con- 
verged solutions without involving an expensive time consuming trial 
and error procedure; 

iii) • For problems with highly nonlinear kinematic behavior, little is under- 
stood of the process of thermomechanical interaction; and lastly, 

iv) Thermomechanically induced pre and postbuckling behavior exhibits in- 
definite stiffness characteristics [2]; such behavior precludes the 
use of the classical form of the incremental Newton Raphson (INR) 
scheme which is restricted to problems with a given definiteness [2,3]. 

Since numerous thermomechanical problems fall into the foregoing cate- 
gories, this paper will consider the development of a solution strategy which 
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bypasses the difficulties denoted by items i) - iv) noted earlier. Specific- 
ally, a constrained type strategy [4-7] will be developed for use with either 
the finite element [1] or difference methodologies. The generality of the 
procedure is such that both pre and postbuckling behavior can be handled 
along with arbitrary kinematic and material nonlinearity. In this context, 
problems exhibiting indefinite stiffness characteristics can be handled. 


GOVERNING EQUATIONS: MECHANICAL 


Assuming the possibility of large deformations, the equations of motion 
complementing the thermal formulation are given by the expression 
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where g 01 - designates the body force vector, <5n< is the Kronecker delta, Sj-j 
the second Plola Kirchhoff stress tensor, ui the deflection vector and aj the 
Lagrangian coordinates. For the current purposes, the Lagrangian strain 
measure Lij is employed in conjunction with s i j namely 
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In terms of the S-fj and Lij measures, the thermoelastic-plastic be- 
havior is handled in terms of the usual yield surface flow rule assumption. 

The creep effects will be treated in terms of strain hardening concepts 
wherein variations in creep rate depend on the existing strain rate. From a 
computational point of view, the overall thermoelastic-plastic-creep behavior' 
is solved via incremental type flow rules. Under the condition of large de- 
formation moderate strain behavior and the usual flow rule assumption, the 
following incremental type constitutive relation is adopted, that is [8] 


AS = [O ep ](AL - AL C - AL t ) (2.3) 

where [Dep] Is the elastic-plastic material stiffness and aL, aLq and ALj 
are increments in Lagrangian creep and thermal strain. For'the'current'work, 
ALc is expressed in terms of mechanical equations of state. In particular, 
it takes the form 


AL C - At y S d 

where S<j is the deviatoric stress and 
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such that o<j and jp- are respectively the equivalent stress and creep strain 

rates. Lastly, the Increment in thermal strain appearing in (2.3) is de- 
fined by 



( 2 . 6 ) 
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where g Is the thermal expansion coefficient matrix and aT is the temper- 
ature increment. Note, based on the thermal fields generated earlier, it 
follows that the various coefficients are temperature dependent. 

In the context of (2.3) it follows that depending on the load step, 
the current stress state is given by the expression 


S = EAS 

where the matrix S takes the form 


(2.7) 
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Noting the linear structure of a§, we see that the incremental scheme 
enables the following segregation of contributing components namely 
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As will be seen later, such a partitioning of the stress state will enable 
the establishment of an improved control of successive iterates during the 
incrementation process. 

FE FORMULATION/SOLUTION ALGORITHM: MECHANICAL 

Following the thermal formulation, we shall employ a displacement 
type procedure to develop the requisite mechanical FE expressions. In this 
context, the deflection field is approximated by 

y = CNylY (3.1) 


where 
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such that [N..] is the displacement type shape function while Y is the nodal 
deflection vector. For consistencies sake, the same order polynomial is used 
for both the thermal and mechanical phases. Based on (3.1) and the virtual 
work principle, the following FE expression can be developed [1] 


M + { dv = F ext 
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Note F noc | a i represents the externally applied nodal 

loads. 

Since dynamic postbuckling problems will be the subject of another 
paper, for the current purposes, we shall consider quasi -static thermo- 
mechanical problems. In this context, (3.3) reduces to the form 

/ K]' S dv - F ext 
v o 

(3.8) 

To simplify the development of the requisite solution algorithm, the 
partitioned form of S will be used to recast (3.8) into a more tractable 
form. Before doing so, we note that due to their analytical form, the 
creep and thermal partitions of S can be lumped with Fext to yield a 
pseudo applied force field namely 

l - -ext * { Kl' ‘!epc + S .epT> dv 
v o 

(3.9) 

Hence (3.8) reduces to the form 
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Since S e pc and S e pT are time dependent terms, the solution to (3.10) 
requires the introduction of a time stepping algorithm to generate the 
requisite solution. This is achieved by expanding (3.10) in truncated 
Taylor series. To start, Y(t+At) is expanded to yield 

Y(t+At) = Y(t) + AY (3.11 ) 
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Substituting (3.11) into (3.10) and truncating higher order terms yields 
the expression 



where 



(3.13) 


such that [S(t)] is the prestress matrix at time t. Based on the defi- 
nition of pseudo force, Eq. (3.9), it follows that 
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Now in terms of (3.10), (3.12) and (3.14), we obtain the following 
time stepping Newton Raphson type algorithm, that is 
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Based on the use of such a relation, successive time steps lead to the 
following thermomechanical history namely 
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As noted earlier, the NR base of (3.15) suffers from several short- 
comings. The more important of these are: 

1. Cannot handle turning points (buckling); 

2. No direct control on successive iterations; and, 

3. Difficult to ascertain zones of convergence as solution proceeds. 


Such drawbacks will be circumvented through the use of constraints in 
the manner of Padovan and Arechaga [6]. Specifically the load increments 
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associated with successive time steps will be constrained. Such a process 
leads to nonuniform time stepping. From the nature of (3.15) it follows 
that constraints must be imposed on increments in the pseudo load F. For 
the current purposes the hyper-elliptic constraint surface (HECS) of 
Padovan, Tovichakchaikul and Arechaga [7] will be employed to control succes- 
sive iterations of a given time step. Such a process is illustrated in 
Fig. 1. The development of the requisite constraint algorithm for the given 
problem requires several main steps, namely: 

i) Establish form of INR extrapolation for a given iteration; 

ii) Establish shape and size of HECS; 

iii) Determine intersection of HECS and INR extrapolation; 

iv) Establish iterative/time stepping aspects of solution algorithm; and, 

v) Establish information required for next time step. 


To start the development, it follows from Fig. 1 that the hyperline 
defining the INR extrapolation takes the form 


On solving for y we obtain 

y = y B + C K UB 3 

such that 

!b = !b " !a 

*B " ?B - Ia 

The HECS appearing in Fig. 1 is given by the following normed poly- 
nomial expression 
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llfll 2 + u A ||y|| 2 = ||f c M 2 
such that 
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The parameter u A appearing in (3.21) regulates the aspect ratio (abscissa/ 
ordinate) of the HECS. 

The intersection of the HECS and INR extrapolation occurs at point I 
as defined in Fig. 1.. Specifically the coordinate? of position I are 
given by 


*i s Ii ’ Ia 

!i = x * (!c " !a } 


(3.22) 
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such that a 1 is a single parameter constraint on the allowable load step 
size and hence the interval in time utilized. Based on (3.22) and (3.23), 
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it follows that (3.17) and (3.20) yield the expressions 


Ki = y .a + + 

(3.24) 
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(3.25) 

In terms of (3.24), (3.25) takes the form 


lu’fcll 2 + u A ||y B ♦ [K u8 r 1 (x I f c - f e )||2 ■ I K c 1 1 2 

(3.26) 

Expanding (3.26) and collecting like terms in A* yields 
polynomial identity namely 

the following 

U 1 ) 2 ^! + 2A I a 2l + « 3 i = 0 
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Solving (3.27) for A*, we obtain 


aI = 5“ { ' a 2l i /(« 2 i) 2 - «u<* 3l > (3.31) 

Based on (3.22), (3.24) and (3.31), Yl the nodal deflection associated with 
the Ith intersection of the INR extrapolation and the HECS takes the form 


ll ■ ?A + *B * ["us]' 1 ( ^7 <-« 2 I i|/(« 2 i) 2 - et i I °3i ) !c - !b> (3 - 32) 

To establish the requisite time stepping aspects of the solution algo- 
rithm, the following variables must be redefined in terms of incrementation 
namely Yj\, Yb, yj-yg, F/\, Fg, Fq and [K(jb]. Letting z denote the time step 
number and i the iteration count, it follows that positions A, B and I in 
Fig. 1 designate the location of the 0th, i th and (i+i)th iterations. In 
this context, it follows that 
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The time associated with point A is the summation of the % preceeding con- 
strained time steps. Hence, since the interval utilized over a given load 
step is xfat, it follows that the time at the end of the «,th step is given 
by 

i „ 

3.36) 


t! s i x7at 
i k=l * 


such that is the finally converged value of the constraint for the k th 
step. K 


Employing the foregoing nomenclature, it follows that 
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dv (3.39) 


The various stress components appearing in (3.37) - (3.39) take the form 
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that the various increments are given by 


4 ?ep<li + 1> 

■ CO ep (vi +1 )]AL(V^) 

(3.46) 

* 

- !t+l> 

(3.47) 

4 5epC^k-l 

+ A k at) ■ - Co ep (v k )]al. c (t k _j ,k k at) 

(3.48) 


160 



“ept^k-l + X k i -‘> * - [t>ep(V k )]«(T(tJ))(T(tVx^t) - T(t X )) (3.49) 

4 !lpC (t X> * ' CV!l + l>Hc (t i- 4t > < 3 ' 50 > 

4 ?ept (t XW ’ - CDepql +1 )]«0-(t x )}(T(t x +i1:) - T(t x )) (3.51) 

To check the convergence of the foregoing algorithm, several tests are 
employed. These include: 

i) Definiteness check: 

( a 2 i) 2 " a ii a 3i > 0 (3.52) 


11) Pseudo force norm check: 
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The proceeding tests are applied at different phases of the iteration 
process. Test i) is used to resize the HECS by self-adaptively readjusting 
u/\ the aspect ratio so as to guarantee an intersection with the INR ex- 
trapolation and thus ensure a convergent solution [6]. Test ii) is employed 
to monitor the monotonicity of successive load excursions. Lastly tests 
ill) and 1v) are used to quantify when adequate convergence has been achieved. 


Once convergence is obtained for a given time step, the overall solu- 
tion algorithm must be prepared for the next interval. This requires that 
the various field variables are properly updated. Specifically this in- 


cludes such terms as Y, Sep. SepC. S e pT and t. In this context, if we let 
I £+1 designate the number of iterations required to yield convergence of 
the U+l)£h time step, then Y at the outset of the U+ 2 )th is given by the 
expression T 
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Note as the iteration process converges, the constraints xj+i represent a 
sequence which approaches the limit (cluster) point X£+i namely 

.... *!., (3-57) 
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In this context, the time at the start of the (*+ l) th step is given by 

‘ui ■ ‘i + X I+ i 4t <3 - 58) 

Now, based on (3.56) and (3.58), it follows that the various stress 
partitions take the form 
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Note for the present purposes, to enhance the speed of calculation 
of the stiffness inverse, the BFGS [3,7] scheme is employed. This 
approach was chosen over the straight updating scheme which is particu- 
larly expensive when several iterations are involved. Such situations 
typically occur in the vicinity of buckling points. 

As was noted earlier, if the INR type scheme is employed to solve 
the thermomechanical problem, uniform time stepping In the thermal phase 
of calculations also leads to equal time intervals for the mechanical 
stage. In contrast, the use of constraints in the INR methodology yields 
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unequal time stepping requirements for the mechanical phase. Namely, 
the following type thermomechanical history Is obtained, that is 
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where here the sequence 0, A, 7 a t,...t* + xJ+iAt is typically nonuniform. 
Because of this, the temperature data required to generate the thermal 
strains and material properties are Interpolated from the uniformly gener- 
ated data. 


BENCHMARKING 

In the preceeding sections, a specialized HECS constrained BFGS up- 
dated INR time stepping strategy has been developed. The methodology 
enables the static solution of pre and postbuckling thermomechanical prob- 
lems. In order to thoroughly evaluate the procedure, several highly non- 
linear benchmark problems were undertaken. The main thrust of this work 
was to ascertain the capability of the constraint methodology to deal 
with thermomechanical problems involving: 

a) Large deformation kinematics including the possibility of pre 
and postbuckling behavior; 

b) Thermoelastlc-plastlc-creep material behavior; 

c) Temperature dependent thermomechanical material properties; 
as well as, 

d) Time dependent thermomechanical loads with varying combinations/ 
Interactions between the thermal and mechanical components. 

This was achieved by programming the solution scheme into ADINA [9] 
and its complementing thermal code ADINAT [10]. Such an approach enabled 
benchmarking over a wide variety of geometric configurations, material types 
and boundary conditions. For the present purposes, the demonstrational 
benchmarking consists of calculating the pre and postbuckling response of 
an arch to various types of thermomechanical loading fields. 

For demonstration purposes. Fig. 2 illustrates the geometry of the 
centrally loaded arch used for the benchmarking. The creep law employed 
Is given by the expression 


As seen from Fig. 2, eight noded plane stress isoparametric elements are 
used to generate the FE simulation. 
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To demonstrate the numerical efficiency and stability of the improved 
constrained MINR time stepping scheme, the thermoelastic-plastic-creep pre- 
postbuckl ing problem depicted in Fig. 3 is considered. As can be seen from 
this figure, the problem is driven into the postbuckling range of behavior 
by the time dependent growth of creep. Overall, the creep generated re- 
shaping initiates a redistribution in the internal loads hence causing a 
change in load carrying capacity. Due to the nature of redistribution, 
plasticity is initiated in the later stages of postbuckling. Noting Fig. 3, 
tcritical marks the tlme at whlch the P re t0 postbuckling transition occurs. 
This time zone is marked- by changes in the definiteness of the structural 
stiffness. Table 1 illustrates the numerical efficiency/stability of the 
BFGS updated constrained scheme in capturing such behavior. In the case of 
At = .8 hours, Table 1, the improved algorithm yielded 210% reduction in 
computer time over the constrained MINR scheme. Note the classical un- 
constrained INR scheme completely fails in such zones of behavior for any 
choice of At. As the time step is increased further, unless some inter- 
mediate updating is employed, even the constrained MINR approach fails. 

This is in contrast to the BFGS updated scheme which shows significantly 
enhanced convergence, efficiency and stability characteristics. 

As a more severe test of the scheme, we shall consider the case of 
cyclical creep loading problems wherein buckling occurs after several cycles. 
Figure 4 illustrates the load deflection behavior of the arch under a 
cyclically applied external load. As can be seen, as the load is cycled 
the accumulated creep over the various cycles progressively reduces the 
buckling limit of the arch. In essence, after several cycles the arch 
behaves as a structure with shape imperfections. Such reductions in 
load carrying capacity are illustrated in Fig. 5. Specifically, this 
figure depicts successive families of load-deflection curves which il- 
lustrate the decrease of buckling strength with time. Note, due to the 
efficiency and stability of the improved constrained MINR time stepping 
scheme, problems involving variable/cyclic loading environments can be 
handled mord effectively. 

The last example considered consists of the thermally induced buck- 
ling of the bimetallic arch depicted in Fig. 6. Noting Fig. 7, as the 
arches temperature is raised, a critical value is reached wherein exces- 
sive deflections occur with no essential raise in T. Such behavior con- 
stitutes the thermal equivalent of buckling. This follows from the fact 
that the structural stiffness is indefinite during the event. 

SUMMARY 

As noted earlier, the main thrust of this work has been to develop an 
improved solution procedure for elastic-plastic creep pre-postbuckl Ing 
problems. Also of major importance Is the maintenance of maximum algo- 
rithmic compatibility with currently available general purpose codes such 
as ADINA, ANSYS, MARC, NASTRAN, etc. As can be seen from the proceeding 
benchmarking, the improved constrained scheme developed herein significantly 
enhances the numerical operating characteristics of MINR type algorithms. 

It should be further noted that due to the manner of formulation, the over- 
all procedure can be encoded into most general purpose codes with little 
rearchitecturing of the programming. 
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Fig.l HECS constrained INR Iterations 
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